/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2024 AUTHOR,AFFILIATION
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    getNearPatchVar

Description

// 1. 获取模型中所有边界名
// 2. 得到靠近具体边界名那一层网格的值（scalar, vector）
// 3. Calculate time averaged shear stress

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<class Type>
Foam::tmp<Foam::Field<Type>> patchInternalField(const UList<Type>& f_in,const UList<Type>& f_bc, fvMesh& mesh,const word& patchName) 
{
    const fvPatch& patch = mesh.boundary()[patchName];

    // Store value on cell near the patch
    tmp<Field<Type>> tpif(new Field<Type>(patch.size()));
    Field<Type>& pif = tpif.ref();
    
    // Store value on face of patch
    tmp<Field<Type>> tpff(new Field<Type>(patch.size()));
    Field<Type>& pff = tpff.ref();

    const labelUList& faceCells = patch.faceCells();
 

    std::ofstream writeDuDy(patchName);
    forAll(pif, faceI)
    {
        pff[faceI] = f_bc[faceI];
        pif[faceI] = f_in[faceCells[faceI]];
        pif[faceI] = pff[faceI] - pif[faceI]; // \delta u

        Info << "At near face cell: "      << mesh.C()[faceI] << endl;
        Info << "Delta u is: "             << pif[faceI][0] << endl;
        Info << "y at face: "              << patch.Cf()[faceI][1] << endl;
        Info << "y at cell near the face:" << mesh.C()[faceCells[faceI]][1] <<endl;
        Info << "Delta y is: "             << patch.Cf()[faceI][1] - mesh.C()[faceCells[faceI]][1] << endl;
        Info << "Delta u / Delta y is: "   << pif[faceI][0] / (patch.Cf()[faceI][1] - mesh.C()[faceCells[faceI]][1] + SMALL) << endl;
        Info << "---" << endl;

        writeDuDy << "at near face cell: ("
                  << mesh.C()[faceI][0] <<", "
                  << mesh.C()[faceI][1] <<", "
                  << mesh.C()[faceI][2] <<") "
                  << "got dudy: ";
        writeDuDy << (pif[faceI][0] / (patch.Cf()[faceI][1] - mesh.C()[faceCells[faceI]][1] + SMALL)) << "\n";
    }
    writeDuDy.close();
    return tpif;
}


int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"

    // 1. 获取模型中所有边界名
    const wordList patchNames = mesh.boundaryMesh().names();
    Info << "The patch names are: " << patchNames << endl;

    // const DimensionedField<scalar, volMesh>& p_internal = p.internalField();
    const DimensionedField<vector, volMesh>& U_internal = U.internalField();
    // const DimensionedField<scalar, fvsPatchField, surfaceMesh>& p_boudnary = p.boundaryField();
    // const DimensionedField<vector, fvsPatchField, surfaceMesh>& U_boundary = U.boundaryField();

    // 2. 得到靠近具体边界名那一层网格的值（scalar, vector）
// Foam::tmp<Foam::Field<Type>> patchInternalField(const UList<Type>& f, fvMesh& mesh,const word& patchName) 
    forAll(patchNames, patchNameI)
    {
        // const fvPatchScalarField& p_boudnary = p.boundaryField()[patchNameI];
        const fvPatchVectorField& U_boundary = U.boundaryField()[patchNameI];
        Info << "------ For patch name: "<< patchNames[patchNameI] << endl;

        // tmp<Field<scalar>> p_nearPatch = patchInternalField(p_internal,p_boudnary, mesh, patchNames[patchNameI]);
        tmp<Field<vector>> U_nearPatch = patchInternalField(U_internal,U_boundary, mesh, patchNames[patchNameI]);

        // forAll(U_nearPatch(), cellI)
        // {
        //     // Info << cellI << ", The value of p near patch " 
        //     //      << patchNames[patchNameI] << " is "
        //     //      << p_nearPatch()[cellI] << endl;
                 
        //     Info << cellI << ", The value of U near patch "
        //          << patchNames[patchNameI]  << " is "
        //          << "(" 
        //          << U_nearPatch()[cellI][0] <<", " 
        //          << U_nearPatch()[cellI][1] <<", "
        //          << U_nearPatch()[cellI][2] << ")" << endl;
        
        // }
        Info << "------" << endl;
    }
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< nl;
    runTime.printExecutionTime(Info);

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //
